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Application  of  Front  Tracking  to  Two  Dimensional  Carved 

Detonation  Fronts 

Bruce  Bukiet  ^^ 

Couxant  Institute  of  Mathematical  Sciences 

New  York  University 

New  York,  N.Y.  10012 


ABSTRACT 

The  method  of  front  tracking  is  applied  to  cylindrically  symmetric  detona- 
tions. This  two  dimensional  method  is  compared  to  a  one  dimensional  ran- 
dom choice  computation  in  the  radial  direction  which  explicitly  utilizes  the 
s>'mmetry  of  the  problem.  The  model  of  detonations  is  tfie  Chapman- Jouguet 
thin  flame  model.  The  random  choice  computation  gives  a  solution  of  high 
accuracy  of  these  model  equations  and  is  thus  taken  as  exact.  We  present 
results  demonstrating  the  high  quality  of  the  front  tracking  solution,  even  on 
coarse  grids.  The  solution  at  the  front  and  in  the  interior  converges  linearly 
to  the  exact  solution.  The  computation  in  the  interior  used  a  second  order 
method  (Lax-Wendroff)  and  the  linear  convergence  appears  to  result  from 
the  fact  that  the  solution  in  the  interior  (i.e.  excluding  discontinuities)  is  C^ 
but  not  C^.  This  is  due  to  "comers"  in  the  solution,  located  at  the  edges  of  a 
rarefaction  wave. 
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Application  of  Front  Tracking  to  Two  Dimensional  Curved 

Detonation  Fronts 

Bruce  Bukiet  ^^^ 

Courant  Institute  of  Mathematical  Sciences 

New  York  University 

New  York,  N.Y.  10012 


1.  IntrodactloD.  The  application  of  the  method  of  front  tracking  to  two  dimensional 
fluid  flows  with  curved  detonation  fronts  is  examined.  We  study  as  test  problems  two  cases 
of  cylindrically  syramietric  expanding  detonation  waves,  namely  a  strong  detonation  and  a  CI 
detonation.  These  problems  can  be  solved  accurately  using  a  one  dimensional  random  choice 
computation  with  operator  splitting.  We  find  that  the  solutions  to  these  problems  using  front 
tracking  converge  linearly  to  the  random  choice  solutions  and  give  satisfactory  answers  on 
moderately  coarse  grids.  We  thus  conclude  that  front  tracking  is  a  valuable  method  for  com- 
puting fluid  flows  with  curved  detonation  fronts. 

2.  Elementary  Theory.  In  this  section,  we  describe  the  elementary  theory  of  one  dimen- 
sional (planar)  Chapman- Jouguet  (CJ)  detonation  waves.  The  equations  of  gas  dynamics  in 
conservation  form  are 

(2.1)  Wt  -I-  Vf(w)  =  0 

where 


w  =    m      and  f(w)  = 


m 
P 


.«_+/. 


^P 


Here  p  is  the  density  of  the  gas,  m  is  the  momentum  density,  and  P   is  the  pressure.  The 
energy  per  unit  volume,  e,  may  be  written  as 
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•2' 
where  v  is  the  velocity  and  c  is  the  specific  internal  energy, 

The  gas  is  polytropic  ('y>l)  and  may  undergo  a  chemical  reaction  in  which  energy  q  is 
released.  We  assume  that  the  chemical  reaction  takes  place  instantaneously  and  thus  that  the 
reaction  zone  is  infinitely  thin.  Suppose  there  is  a  detonation  wave  connecting  constant  states 
labelled  by  subscripts  0  and  1.  From  the  conservation  of  mass,  we  have 

(2.2)  P:(vi  -U)  =  Po(vo  -  U), 

where  U  is  the  speed  of  the  detonation  wave.  From  the  conservation  of  momentum,  we  have 

(2.3)  Pi(vi  -  Uf  +  Fi  =  Po(vo  -  Uf  +  Po  • 
These  relations  lead  to 

P,  -Pn 


M^  = 


where  t  =  —  is  the  specific  volume,  and  M  is  the  the  mass  flux  across  the  detonation  wave. 
P 

For  a  right  wave,  A/^  =  -P]i{v^-Uj^.    For  a  left  wave,  A/^  =  Piiyi-Ui).    These  equations 
are  the  same  as  those  for  shock  waves. 

Equations  (2.2)  and  (2.3)  and  conservation  of  energy,  expressed  by 

(2-4)  (vi  -  U]  [e,  +  i^  J  =   (vo  -  U]  [e,  +  P,] 

yield 


called  the  Hugoniot  relation. 

Suppose  the  gas  in  state  0  is  unbumed  and  that  the  gas  in  state  1  is  burned.  The 
Hugoniot  curve  describes  in  the  (t,  /')-plane  the  states  to  which  the  gas  in  state  0  may  be 
connected  by  a  single  combustion  wave  (See  Fig.  [1]).    The  portion  of  the  curve  bordering 
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the  shaded  region  corresponds  to  non-physical  solutions,  where  M^  <  0.  The  portion  of  the 
curve  to  the  right  of  the  shaded  region  corresponds  to  deflagrations.  The  portion  to  the  left 
corresponds  to  detonations. 

A  line  drawn  from  (tq,  Pq)  that  is  tangent  to  the  detonation  portion  of  the  Hugoniot 
curve  divides  it  into  two  parts.  The  upper  part  corresponds  to  strong  detonations  and  the 
lower  part  to  weak  detonations.  The  point  of  tangency  is  called  the  Chapman- Jouguet  (CJ) 
point.  At  a  CJ  detonation,  the  detonation  wave  moves  at  sound  speed  with  respect  to  the  gas 
behind  it  [4,p.212].   Thus, 

/•,  -  Pn 


(2.6) 


behind  a  CJ  detonation,  where  c^  =  •^^. 

P 

Combining  (2.5)  and  (2.6),  we  get,  after  some  algebra, 
(2.7)  Pc/  + 


2   ,    ,2(<?i-<7o)(7i-l)        2iy-l)P,] 


Pa 


To  7o- 1 

2(y-l)(q-go)Po    ,     (-yi-l)(7o+l) 


=  0 


(7i+1)to  (■yi+i)(7o-i) 

which  is  solved  by  the  quadratic  formula.  Thus,  the  CJ  state  is  completely  determined  by  the 
state  ahead  of  the  detonation.  The  larger  of  the  two  solutions  of  equation  (2.7)  is  the  CJ 
detonation  pressure.  The  smaller  value  corresponds  to  a  CJ  deflagration.  We  get  the  values 
of  pcy  and  t^j  from  equation  (2.6).  Also,  from  (2.6),  M  =  PcjCcj  =  Icj^aPcj-  We  use 
equations  (2.2)  and  the  mass  flux  equations  to  get  v^y  and  Ucj. 

3.  Solution  of  One  Dimensional  Riemann  Problems  with  Combostion.  Consider  equa- 
tion (2.1)  with  the  initial  data 

Sl  =  (Pl.^lPl)  for  j:  <  0 

Sr  =  (pR .«« A)  ^°^  ^  >  0. 

The  solution  at  a  later  time  will  consist,  in  general,  of  a  right  wave,  a  left  wave  and  a 
contact.  We  allow  the  right  (or  left)  wave  to  be  a  shock,  rarefaction,  strong  detonation,  CJ 
detonation  or  CJ  detonation  followed  by  a  rarefaction.  Across  the  contact,  the  pressure  and 
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velocity  are  continuous,  though  the  density  may  jump.  In  order  to  include  at  least  some 
aspects  of  quenching  in  the  model,  we  introduce  a  critical  temperature  (T^)  and  allow 
combustion  to  occur  only  if  the  midstate  temperature  is  above  T^. 

We  write  the  equations  for  the  states  to  which  we  may  connect  the  left  and  right  states. 
Since  the  left  and  right  states  connect  to  states  with  the  same  velocity  and  pressure,  the  equa- 
tions are  expressed  as  velocity  as  a  function  of  pressure.  Then  the  Riemann  problem  is 
solved  by  finding  the  intersection  of  these  curves.  To  find  the  equations  for  shocks  or  deto- 
nations, eliminate  p^  and  U  from  eqs.  (2.2)  (2.3)  and  (2.4),  letting  9i=9o  f°f  shocks.  The 
equations  for  r£u:efactions  are  found  using  the  Riemann  invariants 

Vr  Ct  V.  C, 

-r-  +  r  =  -^  "^  T      for  a  left  rarefaction 

2         -/L-l         2        -y.-l 

Vj>  Cj,  V.  c. 

— r  =  -;; r      for  a  right  rarefaction 

2         7JJ-1         2        7.-1 

as  well  as  the  isentropic  law 

PiPi"^'-  =  P-p."***      on  the  left 

PjtPn~^'  =  F.p."^*      on  the  right, 

where  ♦  denotes  the  state  on  the  inside  of  the  wave.  For  a  CJ  detonation  followed  by  a  rare- 
faction, we  find  the  CJ  state  and  use  the  equations  for  rarefactions,  substituting  CJ  for  the 
subscript  L  or  R. 

We  use  a  Newton  method  iteration  scheme  [1]  to  find  the  intersection  point  and  thus 
solve  the  Riemann  problem.  The  algorithm  initially  finds  a  solution  with  only  hydrodynamic 
waves.  If  the  hydrodynamic  midstate  temperature  (T=P/^)  is  above  some  prescribed  critical 
temperature,  i.e.  the  temperature  at  which  the  unbumed  gas  bums,  and  the  gas  is  unbumed, 
then  the  hydrodynamic  solution  is  replaced  by  a  solution  utilizing  combustion  waves.  The 
random  choice  method  [5]  has  been  applied  to  reactive  gas  dynamics  in  one  dimension  [3], [6] 
using  the  above  solution  of  this  Riemann  problem. 
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4.  The  One  Dimensional  Compotitlon  with  Cylindrical  Symmetry.  With  cylindrical 
symmetry,   the  equations  of  gas  dynamics  (2.1)  become 


Wj  +  f(w)^  =  -S  , 


where 


S  = 


m 

r 

pr 

—{e+P) 
IP*" 


and  r  is  the  radius. 

Sod  [7]  developed  a  method  for  solving  this  system  using  the  random  choice  method  and 
operator  splitting.  At  the  start  of  a  time  step,  we  solve  the  homogeneous  system 

(4.1)  Wj  +  f(w)^  =  0 

by  the  random  choice  method.    Next,  we  use  the  solution  of  eq.  (4.1)  as  initial  data  for  the 
system  of  ordinary  differenial  equations 

w.  =  -S. 


For  p  ^  0,  this  can  be  simplified  to 


p 

r 

V 

=  - 

0 

p 

7Fv 

.   ■ 

t 

I    ''    J 

which  is  solved  analytically.   We  use  a  planar  CJ  thin  flame  detonation  as  described  in  53,  in 
the  cylindrical  solution.  This  approximation  is  exact  in  the  limit  of  large  radius. 

5.  The  Tto  Dimensional  Front  Tracking  Method.  We  present  a  brief  summary  of  front 
tracking  [2]  which  models  two  dimensional  fluid  flow.  A  front  consists  of  a  system  of 
curves.  These  curves,  which  are  usually  associated  with  physical  waves  of  discontinuity, 
evolve  with  the  solution. 
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The  portion  of  the  computational  region  away  from  the  front  is  called  the  interior.  In 
the  interior,  we  solve  an  initial/boimdary-value  problem,  where  the  curves  play  the  role  of 
boundaries.  In  our  computations,  an  operator  split  Lax-Wendroff  method  is  used  to  solve 
eq.  (2.1)  in  the  interior. 

The  propagation  of  the  front  requires  both  the  movement  and  updating  of  the  states  on 
the  front.  Eq.  (2.1)  is  written  in  terms  of  its  normal  and  tangential  derivatives  and  operator 
splitting  is  employed.  First,  a  non-local  Riemann  problem  for  the  detonation  theory  of  §3  is 
solved  in  the  normal  direction,  propagating  the  discontinuity.  Then,  the  solution  of  the  nor- 
mal propagation  is  taken  as  initial  data  for  the  tangential  propagation.  A  one  dimensional 
Lax-Wendroff  scheme  is  used  for  the  tangential  equations. 

6.  Results:  Strong  and  CJ  Detonations.  In  this  section,  we  present  the  results  of  com- 
putations for  expanding  circular  (cylindrically  symmetric)  and  elliptical  strong  detonations. 
The  circular  geometry  was  chosen  for  our  test  problem  because  it  is  essentially  a  one  dimen- 
sional problem  and  the  solutions  by  front  tracking  can  therefore  be  compared  to  solutions  by 
the  random  choice  method.  The  first  sequence  of  runs  models  a  strong  detonation.  The 
second  sequence  shows  the  transition  from  strong  detonation  to  CJ  detonation.  The  two 
dimensional  computations  were  done  on  the  following  grids:  5  by  5,  10  by  10,  20  by  20,  40 
by  40  and  80  by  80.  In  these  computations,  we  tracked  the  contact  as  well  as  the  detonation. 
The  one  dimensional  computations  used  1500  points  in  the  radial  direction,  and  are  taken  as 
giving  the  exact  solution,  within  the  modelling  assumptions  of  the  planar  strong  and  CJ  deto- 
nation theory. 

In  Figs.  2a-f  we  present  the  results  of  computations  in  which  the  pressure  inside  the  ini- 
tial circle  is  100  times  the  pressure  outside.  The  density  is  uniform  and  the  gas  releases 
92.65%  of  its  internal  energy  upon  combustion.  The  pressure  and  heat  release  were  chosen 
to  yield  a  strong  detonation  from  the  start  of  the  computation  through  the  time  step  analyzed. 
Figs.  2a-e  refer  to  a  front  tracking  computation  on  a  40  by  40  grid.  Fig.  2a  shows  the  posi- 
tions of  the  contact  (inner  quarter  circle)  and  the  detonation  (outer  quarter  circle)  at  the  start 
of  the  rim.  The  initial  pressure  profile  is  shown  in  Fig.  2b.  Fig.  2c  shows  the  positions  of 
the  contact  and  the  detonation  at  the  time  step  analyzed.  The  corresponding  pressure  profile 
is  shown  in  Fig,  2d.    In  this  profile,  the  solid  curve  is  the  result  of  the  one  dimensional 
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calculation,  while  the  vertical  error  bars  show  the  range  of  pressure  values  in  the  front  track- 
ing solution  at  fixed  radii,  indicating  angular  dependence.  In  Fig.  2e,  we  present  the  detona- 
tion wave  speed  plotted  against  time  for  the  one  dimensional  run  (solid  curve)  along  with  the 
range  of  speeds  of  the  detonation  in  the  two  dimensional  computation  (vertical  bars).  We 
note  that  the  front  tracking  detonation  speed  error  is  less  than  %%.  Fig.  2f  shows  the  conver- 
gence of  the  front  tracking  solution  to  the  random  choice  solution.  Since  the  solution  is  not 
smooth,  the  Lax-Wendroff  method  used  in  the  interior  is  only  linearly  convergent.  For  each 
grid,  we  compute  the  pressure  error,  in  percent,  at  the  front  and  in  the  interior.  These 
errors  are  plotted  against  the  reciprocal  of  the  number  of  grid  points  in  the  z  £md  y  direc- 
tions. Since  we  use  the  same  number  of  grid  points  (N)  in  each  direction  and  the  computa- 
tional region  is  .5  by  .5,  the  abscissa  is  -j^  =  2Ax  =  2\y.   The  percent  error  in  the  integral 

of  pressure  (#'s)  over  the  computational  region  (interior  error),  the  range  of  errors  of  pres- 
sure, in  percent,  (error  bars)  behind  the  front  and  the  percent  error  of  the  average  value  of 
pressure  behind  the  front  (asterisks)  are  shown  in  Fig.  2f.  The  lines  show  the  errors  going 
linearly  to  0  as  Aj:  =  Ay  -  0. 

Similarly,  in  Figs.  3a-d,  we  show  results  for  an  initial  pressure  ratio  inside  to  outside  of 
68.  In  these  runs,  the  parameters  were  chosen  to  yield  a  CJ  detonation  for  approximately 
three-fourths  of  the  elapsed  time  at  the  time  analyzed  in  these  figures.  The  tenfold  smaller 
errors  in  detonation  speed  in  Fig.  3b  result  from  the  smaller  variations  in  the  speeds  in  this 
40  by  40  nin.  In  Figs.  3c  and  3d,  we  present  separately  the  interior  error  and  front  error  for 
pressure. 

Finally,  we  present  an  example  for  which  a  one  dimensional  method  cannot  be  used. 
In  this  example,  we  find  hot  spots  behind  the  front  in  regions  of  small  curvature  and  cold 
spots  in  corresponding  regions  of  large  curvature.  We  consider  an  elliptical  expanding  deto- 
nation with  the  same  initial  states  as  those  in  the  strong  detonation  runs  above.  The  initial 
lengths  of  the  major  and  minor  axes  are  .3  and  .15  for  the  detonation  wave  and  .29  and  .145 
for  the  contact.  Fig.  4  shows  pressure  contours  and  the  waves  just  before  the  detonation 
wave  breaks  through  the  boundary  on  a  30  by  30  grid.  We  see  that  the  pressure  is  higher 
behind  the  flatter  portion  of  the  detonation  wave  than  behind  the  rounder  portion  of  the 


wave. 

7.  Conclusion.  We  have  shown  that  for  modelling  fluid  flows  with  q'lindrically  sym- 
metric detonation  waves,  the  solution  by  the  method  of  front  tracking  converges  linearly  to 
the  exact  (one  dimensional  random  choice)  solution.  We  believe  these  results  validate  front 
tracking  for  the  computation  of  two  dimensional  fluid  flows  with  curved  detonation  fronts. 
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Fig.  1.  A  Hugoniot  curve  showing  the  states  to  which  Ae  state  (Tq  Pq)  may  be  connected 
by  combustion  waves.  The  shaded  region  corresponds  to  rum-physical  solutions.  The 
lines  connecting  (Tq,  Pq)  to  the  CJ  points,  and  0uis  tangent  to  the  Hugoniot  curve,  are 
drawn. 
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Fig.  2.  The  strong  detonation,  a.  The  initial  strong  detonation  wave  (outer  circle)  and 
the  contact  discontinuity  for  a  cylindrically  symmetric  computation.  The  initial  condi- 
tions are  uniform  density,  zero  velocity,  and  a  circular  pressure  discontinuity  at  radius 
.2.  with  ratio  inside  to  outside  cf  100.  The  heat  released  upon  combustion  is  92.65%  cf 
the  internal  energy  of  the  unturned  gas.  The  initial  position  cf  the  contact  is  radius 
.195.  The  initial  state  between  the  waves  is  that  behind  a  planar  detonation  wave  with 
the  above  initial  data. 


-12- 


^ 
S 


RADIUS 


Fig.  2b.  A  plot  cf  pressure  vs.  radius  at  time  =  0. 
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Fig.  2c.   The  detonation  wave  and  the  contact  at  the  time  sup  analyzed  in  Figs.  2d  and 
21.  The  detonation  wave  now  has  radius  approximately  36. 
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Fig.  2d.  A  plot  of  pressure  vs.  radius  corresponding  to  Fig.  2c  is  shown.  The  solid 
curve  shows  the  results  obtained  by  the  one  dimensional  random  choice  computation.  The 
vertical  lines  represent  the  range  of  pressure  values  in  the  two  dimensional  front  tracking 
solution  at  a  fixed  radius  as  the  angle  varies  on  a  40  by  40  grid.  Thus,  the  vertical  lines 
show  the  angular  dependence  in  the  solution. 
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Fig.  2e.  A  plot  cf  detonation  speed  vs.  time  for  the  computation  on  a  40  by  40  grid. 
The  solid  curve  shows  the  speed  cf  the  detonation  wave  In  the  one  dimensional  calcula- 
tion.   The  vertical  lines  represent  the  range  cf  values  of  the  speed  of  the  detonation  in 


the  two  dimensional  calculation.   The  maximum  error  (max 


U 


-)  is  less  than  %%. 
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Fig.  2f.  Convergence  of  the  front  and  interior  schemes.  The  pressure  errors  in  the 
interior  and  at  the  front  are  shown  for  NxN  grids  at  the  time  indicated  by  Fig.  2c.  The 
#  signs  represent  the  interior  error,  where 

'n\P2i-Pl<i\dxdy 

0  0 


Interior  Error  =  700%  X  -^-; 

0  0 
The  front  error  (error  bars)  gives  the  range  of  the  errors  at  the  front,  defined  as 

P2d~Pjtl 


Front  Error  =  700%  x 


IPJ 


where  IP]  is  the  pressure  jump  at  the  front  in  the  one  dimensional  computation  at  the 
same  time.  The  asterisks  represent  the  error  of  the  average  pressure  behind  the  front, 
namely 

P  ~  P 

Front  Error  (average  pressure)  =  700%  X  °^pi ' 
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Fig.  3.  The  CJ  detonation,  a.  A  plot  of  pressure  vs.  radius  for  a  run  similar  to  that 
described  in  Figs.  2a-e  but  for  a  pressure  ratio  cf68  is  shown.  The  solid  curve  and  the 
vertical  lines  represent  the  one  and  two  dimensional  results,  as  explained  in  the  caption 
to  Fig.  2d.  The  initial  conditions  are  uniform  density,  zero  velocity,  and  a  circular 
pressure  discontinuity  at  radius  .2.  The  heat  released  upon  combustion  is  92.65%  of  the 
internal  energy  of  the  unburned  gas.  The  initial  position  of  the  contact  is  radius  .195. 
The  initial  state  between  the  waves  is  that  behind  a  planar  detonation  wave  with  the 
above  initial  data. 
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Fig.  3b.  A  plot  of  detonation  speed  vs.  time  for  the  computation  on  a  40  by  40  grid. 
The  solid  curve  shows  the  speed  of  the  detonation  wave  in  the  one  dimensional  calcula- 
tion.   The  vertical  lines  represent  the  range  of  values  of  the  speed  of  the  detonation  in 

the  two  dimensional  calculation.    The  maximum  error  (max         .7^^"^  >  «  less  than 
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Fig.  3c.  Convergence  of  the  interior  scheme.  The  pressure  errors  in  the  interior  are 
shown  for  NxN  grids  for  the  time  represented  by  Fig.  3a.  The  error  is  defined  as  in  Fig. 
2f. 
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Fig.  3d.  Convergence  of  the  front  scheme.  The  front  error  is  shown  for  NxN  grids  for 
the  time  represented  by  Fig.  3a.  The  front  error  (error  bars)  gives  the  range  of  the  errors 
at  the  front.  The  asterisks  represent  the  error  of  the  average  pressure  behind  the  front. 
These  errors  are  as  defined  in  the  caption  to  Fig.  2f. 
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Fig.  4.   Pressure  contours  are  shown  for  a  computation  of  an  elliptical  expanding  deto- 
nation on  a  30  by  30  grid.  Also  shown  are  the  detonation  wave  (D)  and  the  contact  (C). 
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